Sandpile model on a quenched substrate generated by kinetic self-avoiding trails 
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Kinetic self- avoiding trails are introduced and used to generate a substrate of randomly quenched 
flow vectors. Sandpile model is studied on such a substrate with asymmetric toppling matrices 
where the precise balance between the net outflow of grains from a toppling site and the total inflow 
of grains to the same site when all its neighbors topple once is maintained at all sites. Within 
numerical accuracy this model behaves in the same way as the multiscaling BTW model. 

PACS numbers: 05.65.+b 05.70.Jk, 45.70.Ht 05.45.Df 
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Accurate determination of the critical exponents and 
thereby precise distinction of critical behaviors among 
the different universality classes are always regarded as 
very important tasks in the field of critical phenomena 
since this study helps in understanding as well as identi- 
fication of the crucial factors that determine the critical 
behaviors. This problem however is still open in the phe- 
nomenon of Self-organized criticality (SOC) in spite of 
extensive research over last several years. More precisely, 
in the sandpile model of SOC the question if the two 
very important models namely the deterministic model 
by Bak, Tang and Wiesenfeld (BTW) 1] and the stochas- 
tic Manna sandpile 0] belong to the same universality 
class or not has not been fully settled yet. A number of 
works claimed that they belong to the same universality 
class 0- EL 13l where as a number of other papers |E S 
argue in favor of their universality classes being different. 
However, what was very much lacking till recently is the 
precise identification of a key factor which may control 
the two behaviors. 

In SOC H E3 a system evolves to a critical state by 
a self-organizing dynamics under a constant, slow exter- 
nal drive in the absence of a fine tuning parameter. The 
signature of the critical state is the spontaneous emer- 
gence of long ranged spatio-temporal correlations in the 
stationary state. The concept of SOC, has been used 
to explain non-linear transport processes of physical en- 
tities like mass, energy, stress etc. in phenomena like 
sandpiles 0, Bl , earthquakes , forest fires 0] , biolog- 
ical evolution |l3| etc. The transport manifests itself as 
intermittent activity bursts called avalanches. Sandpile 
models are the prototype models of SOC. In spite of ex- 
tensive efforts BTW model resisted to follow the finite 
size scaling (FSS) ansatz and has been shown recently to 
obey a multiscaling behavior [Til Il5j . On the other hand 
scaling behavior in the Manna stochastic sandpile |2| is 
observed to be well behaved 

0,1113. 

A deterministic sandpile model can be defined suit- 
ably on an arbitrary graph by an integer toppling matrix 
(TM) A For example, on a square lattice of linear 
size L, the number of sand grains at site i is denoted 
by hi. Sand grains are added to the system one by one 
as: hj — > hi + 1. A threshold value Hi of the number of 



grains is associated with every site. A toppling occurs at 
the site i when hi > Hi . After the toppling the system is 
updated using the TM of size L 2 x L 2 as: hj — » hj — A 
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for j = 1 to L , where A„ = Hi > for all i and Ay 
for all i 7^ j. Therefore during a toppling at the site i, 
the number of grains at this site is reduced by Ajj where 
as -Ay number of grains flow out to all sites j, j = 
1 to L 2 . A toppling at one site may make some of its 
neighboring sites unstable, which may trigger topplings 
at further neighborhood, thus creating an avalanche of 
topplings in a cascade. The BTW model is a special case 
of deterministic sandpiles where the TM has a simple 
structure like Hi = 4 and Ay = A^ = — 1 for each bond 
(ij), otherwise zero 10] . In the Manna stochastic sand- 
pile model each grain of the toppling site is transferred to 
a randomly selected neighboring site implying that the 
TM has the annealed randomness and the elements of A 
matrix is constantly updated during the whole course of 
a given avalanche. 

Recently a single sandpile model with quenched ran- 
dom toppling matrices is proposed which captures the 
crucial features of different sandpile models [ifj. In this 
model the elements of the TM are quenched random vari- 
ables, once their values are selected in the beginning, 
they remain unchanged. The dynamics of the sandpile 
is followed with this TM and the data for avalanches are 
averaged over different random realizations of TMs. In 
this model, there can be two possible situations. In the 
'undirected' case the TM is symmetric i.e., Ay = A^ 
where as in the 'directed' case the TM is asymmetric 
i.e., Ay =/= Aji in general. Here, Ay is nonzero only 
for i — j and for each bond of the lattice. It is argued 
that the behavior of undirected model is similar to the 
BTW model where as that of the directed model is sim- 
ilar to the Manna model. The distinction between the 
two models is made even more precise by defining two 
quantities like Hi = — Xy=&i Ay i.e., the total number of 
grains distributed to the neighbors in a single toppling 
and H i = — Aji which is the number of grains re- 
ceived by the site i when its every neighbor j topple for 
once. It has been suggested that the precise balance at 
all sites (except at the boundary sites) 



H = H 



(1) 
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ensures that the model obeys the same multiscaling as 
in the BTW model. For the directed model this precise 
balance is absent in general and the model shows FSS 
with the same exponents as in the Manna sandpile model. 

In this paper we extend the results of this paper 
and make the conclusion even more precise. We claim 
that only the precise balance Hi = H i or the absence 
of it determines if a model sandpile would belong to the 
BTW or Manna universality class, irrespective of the TM 
being symmetric or asymmetric. 

A quenched configuration of random flow vectors which 
corresponds to an asymmetric TM whose elements sat- 
isfy eqn. (1) is generated in the following way. Let the 
neighbors of the site i be denoted by 1, 2, 3 and 4. We 
first observe that if one increases the A value of any one 
of the four outgoing bonds, say (z3) by an amount 6, 
the bond (z3) becomes asymmetric and it increases Hi 
by the same amount. Similarly if we increase the A 
value of an arbitrary incoming bond to the site i, say 
(2i) by (5 again, the bond (2i) also becomes asymmetric 
and H i increases by an amount 5. Therefore as a result 
of both the operations the precise balance of Hi = H i 
is strictly maintained. In general a series of such bond 
asymmetrizations can be done randomly by starting from 
any arbitrary site i , selecting randomly an arbitrary out- 
going bond (if), increasing Ay by <5, going to the site j, 
selecting an arbitrary outgoing bond (jk) (=/= (ji)) and 
increasing Ajk also by the same amount <5, then going to 
the site k and so on. The path obviously cannot visit a 
bond of the lattice more than once and the final point to 
stop must be the starting point. Such a path can inter- 
sect itself but always one of the outgoing bonds which has 
not been asymmetrized yet is selected randomly. Since 
at each site on the path the A values of either a single or 
a double pair of incoming and outgoing bonds have been 
increased by the same amount 5 the balance of Hi = H i 
is maintained at all sites on the path. 

A self-avoiding trail is a random walk which docs not 
visit one bond of the lattice more than once Il7j . A ran- 
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FIG. 2: (a) Scaling of KSAT loop length distribution for 
system sizes L = 513, 1025 and 2049. (b) The mean end-to- 
end distance of KSATs after n steps grows as n 2v . 

dom configuration of self-avoiding trail is generated by 
growing a random walk which terminates when a bond is 
visited more than once. In contrast a kinetic self-avoiding 
trail (KSAT) is executed with a little more intelligence. 
At each site, to make a step, the walker first finds out 
the subset of bonds which has not been visited yet and 
then steps randomly along any one of these bonds with 
equal probability. Such a walk can also terminate only 
when it visits the origin for the third time (Fig.l). A 
similar definition of kinetic growth walk or growing self- 
avoiding walks have been studied in the literature and it 
is argued that very long such walks behave in the same 
way as ordinary self-avoiding walks |isj . 

First we observe that KSATs have very interesting and 
non-trivial statistics. For example the probability distri- 
bution that a KSAT returns to the origin for the first 
time after n steps has a scaling form like: 

D(n) ~ L-Pg(n/L d f) (2) 

where the scaling function Q(x) ~ x~~ 7 as x — > such 
that 7 = (3/df and Q(x) — > decreases to zero very fast 
when x — > 1. We estimated df « 1.905, (3 « 2.237 which 
give 7 w 1.174 (Fig. 2(a)). The cut-off exponent df is 
also recognized as the fractal dimension of the KSATs 
since the number of steps on the walks whose sizes are 
of the order of L varies as L df . One can also measure 
the value of df directly. The mean square end-to-end 
distance (R 2 (n)) of the walker from the origin after n 
steps varies as n 2v , where v — l/df. Simulation of walks 
of lengths up to a million steps on a lattice of size L = 
4097 gives v » 0.530 so that d f w 1.886 (Fig. 2(b)). 
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FIG. 3: (Color online) Detailed structure of an avalanche, dif- 
ferent sites have toppled different number of times: 1 (circle), 
2(square), 3(triangle up), 4(diamond) and 5(triangle down). 



Therefore we conclude a mean value of df ~ 1.895. 

KSATs are therefore used to asymmetrize the TM. We 
start with a TM whose all elements are initially zero cor- 
responding to a periodic L x L lattice. The walker starts 
from an arbitrarily selected site, executes a KSAT which 
finally stops when it comes back to the origin for the first 
time. The A values of every outgoing bond visited from 
each site are then increased by <5 which is selected as a 
random integer number between 1 and 2. A number of 
such KSAT loops are then generated one by one start- 
ing from arbitrarily selected sites and with randomly se- 
lected S values. The process stops only when all bonds are 
asymmetrized at least once. The periodic boundary con- 
dition is then lifted. The TM so generated is asymmetric 
in « 92.5 % bonds but maintains the precise balance of 
Hi = H i strictly at all sites except on the boundary. The 
lattice is now ready to study the sandpile model where 
the threshold height at each site is denoted by Hi. Such 
a system has a large fluctuation of threshold heights and 
their average increases with increasing the system size. 

We studied three aspects of the sandpile model on the 
quenched substrate generated by KSATs which are: (i) 
the inner structure of the avalanches (ii) the avalanche 
statistics and the (iii) wave size distributions. We observe 
very close similarities of our model with BTW model in 
all three aspects as reported below. 

Like any ordinary sandpile model, the dynamics starts 
from an arbitrary stable distribution of sand heights and 
then grains are added to the system one by one. The 
system eventually reaches the stationary state when the 
average height per site fluctuates around a mean value 
but does not grow any further. The size of an avalanche 
is measured by the total number of topplings s. In Fig. 3 
we show the picture of an avalanche which has no holes. 
Different sites topple different number of times but the 
set of sites which toppled same number of times form 




FIG. 4: (a) An attempt for the scaling of the avalanche size 
distribution for the present model for L — 128, 256 and 512. 
(b) Comparison of the moment exponents a(q) vs. q and (c) 
da(q)/dq vs. q for the present model (solid line) and for the 
BTW (dotted line) model. 



a connected zone. The avalanche has an inward hier- 
archical structure, the n-th toppling zone is completely 
surrounded by the (n — l)-th toppling zone, with the ori- 
gin situated within the maximally toppled zone. This is 
similar to avalanche structure in the BTW model 



The finite size scaling behavior of the probability dis- 
tribution Prob(s,L) of avalanche sizes has the following 
general form: 



Prob(s,L) ~ L-^T{ 



L D 



(3) 



where the scaling function J-(x) ~ x~ T in the limit of 
x — > and J-{x) approaches zero very fast when x — ► 1. It 
is now known that BTW model does not follow this FSS 
form but has a multiscaling behavior 0,0] where as the 
Manna model follows this FSS behavior quite accurately 
0. 

For the present model it is observed that the collapse 
does not work for a single set of /j, and D and for all val- 
ues of s and L. This is a similar situation as found in 
the BTW sandpile model and also in the case of undi- 
rected quenched model [l6|. For example in Fig. 4(a) 
we have tried an unsuccessful attempt for a data col- 
lapse as: Prob(s, L)L 2A vs. sL~ 2 for L = 128, 256 and 
512. Evidently the three curves separate out from one 
another beyond s/L 2 ~ 1. Even for smaller s values 
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FIG. 5: This flow chart shows that the precise balance 
Hi — Hi or absence of it determines the universality classes 
of different sandpile models. 



within 1 < s < L 2 their slopes differ slightly but system- 
atically as 1.132, 1.135 and 1.144 for L = 128, 256 and 
512 respectively, very similar to BTW model behavior. 

Further to check that the present model indeed be- 
haves like the multiscaling BTW model the various mo- 
ments of Prob are evaluated 0, H3. IT^ . The q-th mo- 
ment of the avalanche size distribution is defined as 
(s q ) = £s«Prob(s,L). Assuming that FSS holds for the 
whole accessible range of avalanche sizes, it is known that 
(s q ) ~ where a(q) = D{q - r + 1) for q > t - 1 and 

a(q) = for < q < r — 1. Estimates of a(q) are ob- 
tained from the slopes of the plot of \og(s q (L)) with logi 
for the three system sizes mentioned above and for 251 
equally spaced q values ranging from to 5. In Fig. 4(b) 
we plot o~(q) vs. q for the present model and compare 
it with a similar plot for the BTW model calculated for 
the same system sizes, the agreement is found to be very 
good, within 2%. For both models a(q) shows marked 
deviation from linearity. To analyze this non-linearity in 
more detail it is an usual practice to calculate da(q)/dq 
which takes the constant value D for large q had the FSS 
been valid. In contrast, in our present case we see in Fig. 
4(c) that da{q)/dq increases steadily with q for q > 1 
and this plot coincides within the same accuracy with a 
similar plot for the BTW model. 

In a stable configuration if a grain is added to a site 
i with a height Hi it topples and the first wave is the 
set of all toppled sites while site i is prevented from a 
second toppling. If % is still unstable after the first wave, 
the second wave propagates. This process continues until 
site i becomes stable and the avalanche stops [2(1 ■ 

The autocorrelation function of the wave time se- 
ries {si, S2, S3, ■ ■ ■} of successive waves |22| is denned 
as C(t,L) = [{s k +tS k }L " (sfc)|]/[(sfc)L - (sfc}|] where 
the (..) refers quenched disorder averaging. This long 
range autocorrelation is the consequence of the coherent 
and uniform spatial structure of each wave. C(t, L) is 
found to grow steadily with L. It scales as C(t, L) ~ 
t~ Tc G(t/L Dc ) with same exponents as undirected model 
with t c w 0.35 and D c w 1. These exponents should be 



compared to 0.40 and 1.02, respectively, as determined 
for the BTW model p^ . 

To summarize, an asymmetric toppling matrix is gen- 
erated using random kinetic self-avoiding trail loops on 
the square lattice. The TMs generated in this way guar- 
antees the precise balance between the outflow of grains 
during a single toppling at a site and the total number of 
grains flowing into the same site when all its neighbors 
topple for once. A deterministic sandpile model is studied 
on such a quenched random lattice and the statistical be- 
havior of its avalanches are compared with that of BTW 
model in a number of ways, namely, the inner structure 
of the avalanches, the avalanche statistics and the wave 
size distributions. Within numerical accuracy excellent 
agreement is observed in all three categories. We con- 
clude, as displayed in a flow chart in Fig. 5, that it is 
only the local flow balance or absence of it, irrespective of 
it being generated from a symmetric or asymmetric TM, 
determines the universality class of the sandpile model. 
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